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Abstract 

The main result of this paper is the derivation of a new expression for the tracer 
subgrid term in level coordinates S^ to be employed in O— GCM. The novel feature is the 
proper account of the random nature of the density field which strongly affects the 
transformation from isopycnal to level coordinates of the variables of interest, velocity and 
tracer fields, their correlation functions and ultimately the subgrid terms. Such an effect 
was neglected in previous work. In deriving our result we made use of measured properties 
of vertical ocean turbulence (Gargett et al . 1 98 1 ) . The major new results are: 

1) the new subgrid expression is different from that of the heuristic GM model, 

2) u ++ (tracer)=^u + (thickness), where u ++ and u + are the tracer and thickness bolus 
velocities. In previous models, u ++ =u + , 

2) the subgrid for a tracer r is not the same as that for the density p even when one 
accounts for the obvious absence of a diffusion term in the latter. The difference stems from 
a new treatment of the stochastic nature of the density, 

3) the mesoscale diffusivity enters both locally and non-locally, as the integral over all z's 
from the bottom of the ocean to the level z. 
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I. Introduction 

It has long been suggested (Montgomery, 1940) that mixing in the ocean's stably 
stratified regions occurs preferentially along surfaces of constant densities, isopycnals, 
rather than along isopotential surfaces. More recently, Lozier wet al. (1994) have shown 
that averaging temperature and salinity over pressure surfaces may give rise to water 
masses that do not exist in the real ocean. As discussed in the previous paper (IV), the 
general procedure is thus to construct mesoseale subgrid terms in isopycnal coordinates and 
then transform them to level coordinates. The problem is far from simple because of the 
random nature of the relation between the density p and the coordinate z. Previous work did 
not fully account for this feature: de Szoeke and Bennett (1992) recognize that "the 
equations of motion, continuity, thermodynamics and scalar conservations averaged 
macroscopically on isopycnals and expressed in level coordinates are not the same as the 
conventional equations Reynolds-averaged at fixed depth", but ultimately they argued 
that in practical applications such a difference is unimportant since p' may be considered a 
small parameter. We reach the opposite conclusion. We show that the first term O(p') 
must also be taken into account since the neglect of this feature leads to some poorly 
defined effective density p( z). For future reference, we shall call it "effective p{ z) 
approximation". The goal of this paper is to take proper account of the random nature of 
the function p(z). The resulting subgrid is different from that of previous work. 

The paper is organized as follows: in sec. II, we study the relations between the mean 
and the fluctuating components of the two fields p( z) and z (p) that follow from the random 
nature of p( z). In sec. Ill, we obtain the expressions for the mean and fluctuating 
components for the velocity and tracer fields in level coordinates in terms of the same fields 
in isopycnal coordinates. In sec. IV, we consider the effect of the random nature of p( z) on 
the subgrid for tracers in level coordinates and obtain a model independent expression for 
the tracer subgrid term in level coordinates in terms of two functions which must be 
modeled in isopycnal coordinates and which we discuss in sec.V. In sec. VI, we present the 
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final result for to be used in O-GCM. In sec. VII we present as the sum of 
diffusion+advection terms which differs considerably from previous expressions. In sec. VIII 
we discuss the difference between our expressions for the bolus velocities for tracer and 
densities in both coordinates with those of previous authors. 


II. Relations between the random fields p{ z) and z(p). 

Consider the two fields p(x,y,z,t) and z(x,y,p,t) which we split in two parts 
representing the average and the fluctuating components, 


Pi z) = P( z) + P\ z ) 

(la) 

Z {p) = Z (p) + z'(p) 

(lb) 

For sake of simplicity, we shall omit the x,y dependence. For arbitrary 

z and p, we have the 

exact relations 


p=p[z{p)}, z=z[p(z)] 

(2) 

Further, we assume that the fluctuating components p'( z) and z'(p) are 

sufficiently small to 

allow a power expansion. Using (3a, b) we then have (the subscripts denote partial 

derivatives) 


p{ z) = p( z+z') = p{ z) + p z { z)z' + \p iz ( z)z' 2 +... 

(3a) 

» z® =f!z=i 

(3b) 

Averaging, subtracting from the original equation and keeping terms up to 0(z' 2 ), we 

obtain: 


p = m + ft 

U) 

where 


fp = ^(z)z' + fp zz { i)z' 2 

(5) 

and 


p'{ z) + p z (z)z'(p) + ¥= 0 

(6a) 

Sp' = A[^(z)z'] + ip zz (z)A[z' 2 ] 

(6b) 
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where we have defined 


A(xy) = xy-xy (7a) 

and in Eqs.(4-6) we have used the shorthand notation 

z=z (/>), z'=z '(p) (7b) 

An analogous procedure with the variable z (p) leads to the following relations which 
correspond to Eqs.(4-5) 


z = z (p) + 6z 


(8a) 


<5z = z' p (p)p' + ^pp(p)p' 2 ( 8b ) 

The counterparts of relations (6) are: 

z '(p) + z p (p)p' + <fe' = 0 (9a) 

&' = A[z^(p)p'] + iz^(p)]A(p' 2 ) (9b) 

where we have used the notation: 

P=P(z), p=p'( z) (9c) 

The "effective p(z) approximation" used by previous authors neglects both terms 6p and 5z 
in (4) and (8a) which become 

p=p\z{p)], z=z[p(z)] (10a) 

Since in these relations p and z are arbitrary, we can choose 

P=p( z) (10b) 

Then, the second of (10a) implies that 

z=z (p) (10c) 

Relations (10b,c) imply that functions z(x) and p(x) are inverse of each other, a property 
that is fulfilled only within the approximation used above. Relations (10) are implicit in 
the work of de Szoeke and Bennett (1993) and frequently used by several authors: for 
example, Eq.(14) of Gent et al. (1995), Eq.(4) of Treguier et al. (1997) and Eqs.(B.4), 
(B.8) and (B.12) of Smith (1999). 

Consider the effect of retaining terms up to O(p') and O(z'). Since Eqs.(4) and (8a) do 
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not contain such terms, relations (10) remain unchanged while Equations (6a) and (9a) 
become: 

p\ z ) + P z ( z ) z '(p) = 0 ( lla ) 

z '(p) + z p (p)p'( z ) = 0 (Hb) 

which proves that 

P z ( z ) z p (p) = 1 (He) 

This relation becomes an identity when one neglects the random nature of p(z . Relation 
(11c) is frequently used and relations (11 a,b) were employed, for example, by Treguier et 
al. (1997) and Rix and Willebrand (1996). 

If we now account for terms 0(p' 2 ) and 0(z' 2 ), in lieu of relations (10a, b), Eqs.(4) and 
(8a) yield: 

p- p[z{p)} = fp (12a) 

z-z[p(z)] = 6z (12b) 

where 8p and 6z are given by expressions (5) and (8b) with (10b, c). Furthermore, using 
Eqs.(10) and (11), we obtain: 

6p = ~P z fc (12c) 

We also write explicitly the second-order approximation in Eqs.(6a) and (9a). They are: 

p'(z) + P z {z)z'{p) =-Sp' =- p z (z)5z' (13a) 

where 

¥ = A[/^(z)z'(p)] + ip zz (z)A[z'(p) 2 ] (13b) 

6z' = A[z^(p)p'(z)] + ±z pp (p)A[p'(z) 2 ] (13c) 

where (10b, c) are assumed. Relations (12)— (13) have never been considered before and yet 
they are key to obtain the correct transformation from isopycnal to level coordinates. 

III. Velocity and tracer fields in isopycnal and level coordinates. 

Due to the random nature of the functions p( z) and z(p), the transformation of 
random fields like velocity and tracer from isopycnal to level coordinates is far from trivial. 
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To begin with, we split the fields into mean and fluctuating parts; in the case of a 3D 
velocity field v, we write 

level: v(z) = v(z) + v'(z) (14a) 

isopycnal: v(p) = \(p) + v"(/>) (14b) 

An overbar represents average in level coordinates, and a tilde represents average in 
isopycnal coordinates; similarly, the fluctuating components are denoted by a prime (') in 
level coordinates and by a double prime (") in isopycnal coordinates. For z (/?), we use 
Eq.(lb) instead of (14b) since the only meaningful average is at fixed p and thus there is no 
need to distinguish between a bar and a tilde. The exact relation between the fields in the 
two coordinate system is written as 

v(z) = v[p( z)J (15) 

Let us substitute in (15) the decomposition (la) for p( z) and let us expand the right hand 
side in power series in p'( z). For the analysis that follows we retain only the zeroth and 
first powers in p'{ z). We have 

v(z) = v[p(z)] + v p [p(z)]p'(z) (16a) 

Substituting Eq.(14b) and separating mean and fluctuating parts, we obtain 

v(z) = v[/?(z)] + V" p [p{z))p\z) (16b) 

v'(z) = v"[p(z)] + \ p [p(z)]p'{z) + Afv^z)]?'^)} (16c) 

where 

v[p(z)] = v(/))|^ (16d) 

v"[p(z )] = v"(/>)| p= -( z ) (16e) 

In analogy with (15), we have the relation 

v(p) = v[z(p)] (17) 

Using (14b) and expanding in powers of z '(/>), we obtain 

v{p) = v[z(/>)] + ( 18a ) 

v"(p) = v'[z(p)] + vjz(p)]z'(p) + A{v^[z(p)]z'(p)} (18b) 
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Notice further that we limit ourselves to terms 0(/?' ) and 0(f' 2 ), where f 1 is any other 
fluctuating variable. Thus, when we average bilinear terms in the fluctuations, as in 
Eqs.(16b) and (18a), the difference between averages in isopycnal and level coordinates is 
negligible provided z=z(p) or p=p( z) since it entails higher order fluctuations. Thus, 
Eq.(18a) can be rewritten as follows: 


*W)\ = v(z) + v' z (z)z'[p{z)\ (18c) 

However, even within this approximation, we must distinguish between the two 
fluctuations v'(z) and v"(p), Eqs.(16c) and (18b). Analogous relations can be obtained for 
the tracer fields r(z) and r(p)\ they can be obtained by substituting 


v->r 

in Eqs.(16) and (18). 


(18d) 


IV. Tracer field. Level Coordinates. Model Independent Result 

\\ hen dealing with the equation for the mean tracer field t in level coordinates (V is 

3 

the 3D gradient operator and v=u,w), 

d t r + v-V 3 r + S^= 0 (19a) 

one must model the subgrid term (VhV h is the two dimensional gradient) 

S^= v '-V 3 r' = u'-Vr’ + w V (19b) 

Due to the fact that it is physically easier to interpret and model variables in isopycnal 
coordinates, we shall first express in terms of correlation functions in such coordinates 
and then transform the result to level coordinates. The first step is the substitution of 
Eq.(16c) and the analogous relation for r'(z) (obtained by substituting v-*r in 16c), into 
Eq.(19b). Since in isopycnal coordinates we only consider horizontal components of the 
velocity field, that is, u(/?), we need to eliminate the vertical components of v(z). To this 

end, we use the continuity equation for v'(z) which yields 

00 

w '( z ) = V-/ z u'(z*)dz* (20a) 
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(20b) 


Transforming the variables of integration so that z*=z(x), we obtain 

p( z ) 

w'( z )=-V-/ 0 u'[z(x)] z^(x) dx 

Substituting the operator V (=V H ) using the relation (V is the gradient at fixed p) 

S = S p + {Hp)d p (21a) 

we obtain 

w'(z) = - ^ z (z)‘ 1 u'(z)-V^(z) 4- tfw'(z) (21b) 

where 

p( z ) 

Sw'(z) = -/ 0 V x -u'[z(x)J (p z ^(x)]}' 1 dx (21c) 

so that 

3>'(z) = -? p u’-u (21d) 

Let us notice that in Eqs.(20b) and (21c) we consider u' and p to be functions of p whereas 
initially in Eqs.(14a) and (la) they were defined as functions of z. The change of the 
independent variables is carried out via the substitution z=z (p). Analogously, even though 
the field v" is initially defined in (14b) as a function of p, in (16c) we express v" in terms of 
z via the substitution p=p( z). Actually, the detailed specifications shown in Eqs.(16c) and 
(21) are not necessary; they can be easily reconstructed since the choices of the independent 
variables are unique. As we discussed earlier, with accuracy up to the second order in the 
fluctuating fields, the average of the bilinear functions of fluctuating fields does not depend 
on which independent variable one chooses, p or z. Thus, for sake of simplicity, we shall 
frequently omit the arguments of the functions and use a overbar instead of a tilde. With 
this clarifications, we substitute Eqs.(21) into Eq.(19b) and obtain 

s ( = + i^' ( 22 ) 

After substituting Eq.(21c), one observes the following: 

1) in the first term there are only one-point correlation functions whereas in the second 
term there two-point correlations functions, 

2) in the first and second term there are different integration and differentiation processes 
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whose relative importance can be estimated by using characteristic horizontal and vertical 
(more accurately, isopycnal and diapycnal) length scales L and H, as well as the 
characteristic diapycnal depth of the correlation, <5H. The experimental results of Gargett 
et al. (1981) show that <5H is of the order of tens of meters which is considerably smaller 
than H>10 3 meters. In this way, we obtain that the ratio of the second to the first term in 
(22) is of order (5H/H and we can neglect it. Thus, we retain only the first term in (22), 

s t-Vf-v 7 *! (23) 

Eq.(23) can also be derived in a different way. Multiply (21b) by r 1 , average and neglect 
the term {wV in accordance with what we have just discussed. As a result, we obtain 

P z wV + uV -Vp = 0 (24) 

Using (24) in Eq.(19b), we derive (23). We further notice that (24) also implies that, the 
subgrid tracer flux v T 7'=(u r F, w'rj is directed along isopycnal surfaces. An similar 
conclusion for the density flux can be obtained from (21b): 

p w'p' + u V - Vp = 0 (25) 

Eqs.(24)-(25) express one of the basic facts of physical oceanography, namely that mixing 
in the stably stratified ocean occurs mostly along isopycnal surfaces. 

Next, we use (16c) and its analog for the field r, to express u' and r' in Eq.(23) in 
terms of u" and r". This yields: 

= ^z v p*^ u " + V ,)(r " + V )/ ^ 1 ( 26 ) 

which we expand as a power series in p'\ 

S/i = S + S + (27) 

« o 1 

where: 

s 0 = v/rrrfj t 2 **) 

To first order in O(p'), we further have: 

* ** 

S = S + S (28b) 

ill 

s * = PjpFpip) - u"(p)p'/p z ) ( 28c ) 
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S i = p z i p -& p (i>W?Wp1 ( 28d ) 

The last term is due to the correlation between the density and the tracer fields which 
interact indirectly through the velocity field. On that basis, we suggest that to order 0 (//), 
this term may be neglected 

S** » 0 (28e) 

* 

Next, we transform S . It is convenient to introduce the velocity (see Eq.lOa of IV) 

u = - d z (uY/p z ) = (z T u TT ) z (29) 

Since at the bottom z=-H, u" vanishes, we have 

uV/P z = - /_ H u (z')dz' = I (30) 

Using Eqs.(27), (28) and (30), expressed in isopycnal coordinates is: 

s <= s „ + ^-(V ) (31) 

One more step is needed to write (31) explicitly in terms of level coordinates. This is 
accomplished by the substitutions 

d p^P'z d z' V p ^-(Vp)P~z d z ( 32 ) 

Thus, Eq.(31) takes the form 

= Vz + (Vp z ) u -Vp + S* + S q (33a) 

where 

W=VI, S, =I-a z [r z (L^-L r )] (33c) 

L f, = ~ <•■>' L T =-/>-'Vr (33d) 

where T are density and tracer surfaces slopes. The variables u^, I and are given by 

(29), (30) and (28a). One must model two variables 

u , S (34) 

10 

which we consider next. 

V. Interpretation of u^ 

Following Andrews and McIntyre (1976) and Andrews et al. (1987), the transport 
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(35a) 


equation for the mean density p is given by: 

d.p + (u+u + )*Vp + (w+w + )p„ + E = 0 

L Z Z 

where 

u + =-5 z A, w + =V-A, A=(^ z )- 1 Uy (35b) 

E = A- Vp + w'p' (35c) 

Consider also the equation for the density variance in the adiabatic case, 

d.p' 2 + 2u‘p' • Vp + 2w '// p + u^p 12 ) + wd (p l2 ) + d (w'p' 2 )+ V-iTp' 2 = 0 

z z z 

(35d) 

If we neglect terms 0(p' 2 ), only the second and third term in (35d) survive, which leads to 
Eq.(25) which in turn implies that 

E =0 (35e) 

Consequently, Eq.(35a) becomes 

d p + (u+u + )-Vp 4- (w+w"*")p = 0 (35f) 

where u + can now be interpreted as a density bolus velocity in level coordinates. We can 
further notice that using (29), (35b) and (18b) we obtain 

u - u + = d z (u z p ,2 /p 2 ) (36a) 

that is, the difference is of higher order than and u both of which are O(p'). Since 
Eq.(33) and (35f) are also correct up to O(p'), in (29) we can substitute 

u"-*u' (36b) 

which implies that in (33a) we can substitute 

= u + (36c) 

Thus, the density bolus velocity in level coordinates u" 1 " coincides with u^ defined by (29). 

VI. Tracer subgrid in Level Coordinates 

Eq.(33) contains two functions that need modeling, a velocity and a diffusion-like 
term S . In paper IV, we modeled u^ and its form is given by Eq.(26). Before substituting it 
into (33), we must transform it to level coordinates. Using Eq.(32), it becomes 
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u ,= K + (37a) 

where we have taken into account that in the case of coarse resolution O-GCMs, we can 

neglect £ compared to f, as well as the fact that f =0. To derive S , we first recall that the 

z o 

general result (33) is valid to second order in the fluctuating fields and to first order in p' . 
Within this accuracy, we can use Eq.(27) of IV which we compare with Eq.(28a) to 
conclude that: 

S q = R (37b) 

where, within the same approximation, R which is given by Eqs.(8a) and (28) of IV, 
acquires the following form 

R= -^y(yy < 37c ) 

where r is the thickness weighted average in isopycnal coordinates, see Eq.(4a) of IV, which 
consistently with the approximation we have used, can be substituted with r that is, 

R = R = “A V (37d) 

Transforming this expression to level coordinates using Eq.(32), we obtain 

<*«) 

where {p-sdp/dx^) 

= Vw 2 (37f) 

Thus, substituting Eqs.(37a,d,e) into (33a), we can present result in two alternatives forms: 
(p=p, tbt): 

1) first form: 


S/t) = v/T z + ( tJp z ) u -Vp + S* + R 

(38a) 

u = nd z {p z Sp) + Kf _1 Vf, w =V-I 

(38b) 

S * = 1 -WVM 

(38c) 

V- N 

1 

II 

k 

1 

S3 

1 

II 

(38d) 

1 = - /_ R U j( z ') dz ' 

(38e) 


(38f) 
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= ^-(i-O^Yk 2 ( 38 §) 

where we have added the diapycnal component with diffusivity kc which is the vertical 
diffusivity. The first two terms in (38a) cannot be interpreted as an advection, as in the 
Gent and McWilliams model (1990, GM). However, we can rearrange (33a) to exhibit such 
term. To this end, we rewrite (37d) as follows: 

R = R*- wjy p P z l )-V p T = R* + u -S p T- rf^Vf-V r (39a) 

where 

R * = " ~ ( 39b ) 

Using (21a), we obtain: 

2) second form: v^= (u ,w ) 

S^r) = + S* + R* - /cT’Vf^V+L d z )r (40a) 

R* = - (V + L p d z ) • «(V +L p d z )r- ^(^P\Pf^^) (40b) 

The expressions for u^, S^, w^ and I are the same as in Eqs.(38). 

Since the first term in (40a) looks like an advection term whereas the first term in 
(38a) does not, two questions arise: first, is a division of the level subgrid in advective plus 
diffusive parts unique representing a physical reality or is just a formal mathematical 
reformulation?, second, if the answer is positive, is the first term in (40a) the physical 
advection term? We discuss this problem in the next section 

3) the GM model 

If we take a constant diffusivity k as GM did, we obtain from (38) 

S / r ) = S/GM,r) + [r z Iy (L^-Lp] (40c) 

which still does not coincide with the GM model. Only if we furher take r=p, will the last 
term vanish and the two models coincide. 

VII. Diffusion and Advection in Subgrid Modeling 

To discuss the problems just mentioned, we recall that following Taylor's (1915) 
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(41a) 


ideas, Rhines and Holland (1979, Eqs.4-9a) have shown that equation (19a) (r=r) 


Df d / — i — r\ o 

ut ®,(v) = _s < 


can be written quite generally as: 


5 P’— “u®. 


S ( - “ lx-( D ij 

1 J 


(41b) 


where the diffusivity is a second-order tensor the orientation of which is unknown. As 
stressed by Muller and Holloway (1989), Greatbatch (1998), the principal axes may be 
taken to be diagonal in level, isopycnals or any other coordinate systems. For a number of 
years, the tensor was assumed to be diagonal in the (xyz) space giving rise to 
horizontal/vertical diffusivities. Redi (1982) suggested that it is diagonal in the 
isopycnal/diapycnal system and since mesoscale mixing occur along isopycnals, the tensor 
has ever since been taken diagonal in that system and thus one refers to diapycnal and 
isopycnal diffusivities. Suppose that D- j is a symmetric tensor and that we take it equal to 
Redi's Kjj, Eq.(38g): 


D..=K.. 
ij iJ 


(41c) 


where 


K ij = *< V Wfk + V 'Yk 2) = K ?j + <K ij (D) (41d) 

The first two terms represent the isopycnal contribution since K9.p.=0, while the last term 


represents the diapycnal contribution since Kjj(D)p.=p i . The question is: can the velocity 

(41e) 


v ++ = |-K.. 
' *j IJ 


be interpreted as the bolus velocity v| + (u| + ,w| + ), see Eq.(lb) of paper IV? The 
horizontal and vertical components are: 


~£^z 1V/> ) + -> (4 If) 

which have the same functional dependence of the GM model bolus velocity but with the 
difference that both velocities (4 If) are either from low to high or vice-versa, whereas in 
the GM model the horizontal part is from low to high f while the vertical is from high to 
low f. In more mathematical terms, the true bolus velocity is divergence free (by 
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construction) while (41f) is not. Next, we write D- j as the sum of a symmetric part, f - and 
an antisymmetric part a^ (Plumb, 1979): 


Thus, 


D- — K- a- • 
ij ij ij 


+ 4c a ij ~ij ra ij) 
j j J j J 


(42a) 


(42b) 


Contrary to (41e), we now consider d/dx-^ of a^ as a possible candidate for the advective 
velocity v"j~ + 

v T + = %x a ij ( 42c ) 

j 

Since a^-a^, it follows that the last term in (42b) is such that 


(42d) 


~ acj^tj^ji) - 0 (42d) 

and thus it does not contribute to in (41a). By the same token, it follows that v"t + is 
divergence free 


Vr = ° 


(42e) 


and so going back to (41b) we have 


S ^ _v ( 42f ) 

This shows that the symmetric part is a diffusion process while the antisymmetric part 

gives rise to an advection velocity. Thus, if D- j is known somehow, the decomposition of S^> 

into advective and diffusive terms is unique since for any tensor the decomposition in 

symmetric and antisymmetric parts is unique. If, on the contrary, D- ■ is not known but 

B * 

is, as in our case where we have constructed Eqs.(38) and (40), the uniqueness of the 
decomposition of S ^ is predicated on D-- depending only on p but not on r. First, let us use 
(42a) into (41b). We obtain (where commas indicated derivatives): 

S l =- K ij r ,i j - < K ij,i + a ij,i > r ,j (42g) 

By equating the terms in the known which contain the second derivatives of r with the 

first term in (42g), one can construct the tensor To obtain the last two term in (42f), 
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we equate them with the corresponding terms in S^, which contain the first derivative of r. 
This will yield a-. .. In summary, by knowing S f we can uniquely reconstruct from it both 
the tensor as w'dl as a- - • which is nothing but — v as Eq.(42c) shows. Applying this 

jj 1J,1 J 

method using the specific form (38) and/or (40), we obtain the following results: 


where: 


s / r > = s diff (r) + W 7 ) 
S difl( T ) = - Ix^ij 

S adv< r > = * + ,% 

(43a) 

(43b) 

(43c) 

v++ _ & o 

' > 

(43d) 

u ++ = 
i i 

(44a) 

Vu ++ + w ++ =0 
1 

(44b) 

< — . . 

II 

J* 

+ 

(44c) 

A ft/? =0, (a,0= 1,2) 

(44d) 

A = A = - £1 , (a=l,2) 

(44e) 

A = p~^ I • V/9 

33 Z 

(44f) 


Thus, we reach the following conclusions: first , separation (43a) in terms of advective and 
diffusive components is unique; second , the velocity in the first term in (40a) is not the 
physically correct bolus velocity since it is twice as large as the real (advective) bolus 
velocity (u'|' + ,wt + ), Eqs.(43a)-(44a); third , the diffusion tensor is not only the Redi 
Kjj, as assumed thus far, since it contains the additional term A.y Since the velocity u^ 
equal the density bolus velocity, Eq.(36c), it may at first seem surprising that in level 
coordinates 


u ++ (tracer) = £u + (thickness) 


(44g) 


instead of 


u^ + (tracer) = u+( thickness) (44h) 

We now show that (44g) and (44h) represent two different interpretations of the physical 
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content of the density subgrid. If we substitute r-*p in (43) we obtain 

S <M = s diffM + s »dvW 

= ~ Ix^ij $*.)> S adv^ = v t + !f ] (44i > 

Using now (43d) and (44a-f), we obtain: 

S diff^) = S adv^) = *( u r V/? + Vz> ( 44 j) 

S/p) = u-V/> + (44k) 

Eq.(44j) can be interpreted in only one way: density subgrid is the sum of advection plus 
diffusion which contribute an equal amount. On the other hand, (44k) can be interpreted in 
two ways: either as the sum of two equal terms or, as the absence of diffusion and full 
contribution of advection. The first choice implies (44h) while the second choice implies 
(449)- Thus, we conclude that while both (44g) and (44h) are acceptable, the common 
interpretation is that for density there is no diffusion and thus (44g) is to be preferred. 
Stated differently, the standard model that: 

s diffM=°. u+ + =u + 
is untenable since the correct relations are either 

s d iff(^) =0 > u t += ’ u+ 
or: 

S diffM = S adv Wi(0 ' “t +=u+ 

The question still remains as to why for density, the separation of advective from 
diffusive terms is not unique. The answer is that in the proof of uniqueness given above, the 
critical condition was that be independent of r which is no longer true when r->p since 
Eq.(44c-f) show that via k - depends on p. 

VIII. Present and previous models 

In principle, there are four bolus velocities: for density and tracer and in level and 
isopycnal coordinates, that is, 


17 



1 * 

density: u (level), u (iso) (45a) 

, , ** 

tracer: u j” t ’(level), u (iso) (45b) 

In terms of these velocities, the subgrid tracer functions ] are defined as: 

S^= u ++ -Vr + w ++ r + R(level) (46a) 

t- I 1 z 

Sj = u -V^r + R(iso) (46b) 

In all previous models it was assumed that all four bolus velocities are identical: 

** * i , , 

u (iso)=u (iso)=u (level)=u"T (level) (46c) 

whereas we have shown that: 


u (iso) = u (iso) (47a) 

u'|’ + (level) = |u"*"(level) (47b) 

u (iso) = u + (level) + (47c) 

where is given by Eq.(24) of IV. To understand how the second equality in (46c) was 
arrived at, we begin with the equation for the thickness z^=dz /dp whose mean value in 
isopycnal coordinates satisfies the equation: 

v [(S+u X ]=o (48a) 

where 


u = u"z ' / z 

p> p 


(48b) 


is the thickness bolus velocity in isopycnal coordinates. Gent and al. (1995) transformed 
(48a) to level coordinates to obtain 


d.p + (u + u )V -/0 + (w + w )d p = 0 (49a) 

Within the O(p') accuracy, in (49a) one takes p=p, a variable for which we have already 
given the transport equation (35f) which we reproduce here for ease of comparison: 

dj) + (u+u + )-Vp + (w+w + )p z = 0 (49b) 

By comparing (49a, b), one can only conclude that 

~ * - + 
u + u = u + u 

w + w = w + w + (49c) 
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whereas in the GM model 


u = u+ (49d) 

which is the assumption we referred to. Thus, we see that (49d) is due to the lack of 
differentiation between u and u [in our model, substituting Eqs.(35b), (48b) and (11a) into 
(18a), one obtains the first of (49c)]. 

To understand whence the last assumption (46c) in the GM model arises, we recall 
that GM heuristically extended (49a) to the case of a tracer by adding a diffusion term R 
thus obtaining (see also Gent et al ., 1995) 


t^r+u-Vr + w r z + S^=0 
S^(GM) = u* • Vr + w*r z + R 

Eq.(50b) is consistent with (46c) only if ones makes the assumption 


(50a) 

(50b) 

(50c) 


IX. Conclusions. 

The goal of this paper was to improve on the "effective p( z) approximation", 
employed thus far in all considerations, since it neglects the intrinsic random nature of the 
function p( z). We began by considering the transformation of the mean and fluctuating 
parts of the different fields of interest from isopycnal to level coordinates. For the velocity 
field, the results are presented in Eqs.(16) and (18) which are accurate up to the 
second-order in the fluctuating fields. Analogous transformations for the tracer fields can 
be obtained via the substitution v-*r. Using these results, we compute the subgrid term 
that appears in Equation (19a) for the mean tracer t in level coordinates. The main results 
are expressed by Eqs. (38) and (40) and/or (43)-(44)- A considerable simplification of the 
problem was achieved thanks to the experimental results by Gargett et al. (1981) that 
showed that the vertical correlation length scale for ocean turbulence (below the mixed 
layer) is much smaller than the characteristic vertical extent, H. As we have shown, this 
feature leads to Eqs. (24,25) which in turn imply that ocean mixing in stably stratified 
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regions occurs mainly along isopycnals. On this basis we derived our main results which are 
quite different from the GM model, Eq.(50b). 
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